function c = bnet3(dofs,mesh,x,f);
%        c = bnet(dofs,V,T,posT,d,f,varargin);
% This function returns the coefficients in the vector c of the piecewise interpolation
% polynomials of the function f(x,y) over the triangulation [V,T]. That is, c is the
% coefficient vector of continuous spline function of degree d which interpolates f(x,y)
% over the domain points of each triangle in triangulation [V,T].
V=mesh.V;
T=mesh.T;
posT=size(T,1);
d=mesh.d;
c=mesh.c;
mat = vdm21(d); m = (d+1)*(d+2)/2; m2 = m*m*posT;
Indx1 = zeros(m2,1); Indx2 = zeros(m2,1);
S = zeros(m2,1);  [I,J,K]=indices(d);
dim = max(max(dofs))-min(min(dofs))+1;
b = zeros(dim,1); pos = 1; 
[i,j,s] = find(mat); L = length(i);
for tri = 1:posT
    pts = (I*V(T(tri,1),:)+J*V(T(tri,2),:)+K*V(T(tri,3),:))/d;
      Indx1(pos:(pos + L-1)) = dofs(tri,i)';
    Indx2(pos:(pos + L-1)) = dofs(tri,j)';
    S(pos:(pos + L-1)) = s;
    pos = pos + L;      
    c_loc=c(dofs(tri,:));
    x_loc=x(dofs(tri,:));
          b_loc = -(mat*c_loc).^3+3*(mat*x_loc).^2.*(mat*c_loc)+feval(f,pts(:,1),pts(:,2)); 
    %b_loc = feval(f,pts(:,1),pts(:,2)); 
    b(dofs(tri,:)) = b(dofs(tri,:)) + b_loc;
end
A = sparse(Indx1(1:(pos-1)),Indx2(1:(pos-1)),S(1:(pos-1)),dim,dim);
c = A\b;